library(raster)
## Loading required package: sp
library(tidyverse)
## ── Attaching packages ──────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.3 ✓ dplyr 1.0.1
## ✓ tidyr 1.1.1 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0
## ── Conflicts ─────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x tidyr::extract() masks raster::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x dplyr::select() masks raster::select()
library(getlandsat)
library(sf)
## Linking to GEOS 3.8.1, GDAL 3.1.1, PROJ 6.3.1
library(mapview)
## GDAL version >= 3.1.0 | setting mapviewOptions(fgb = TRUE)
library(osmdata)
## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
library(factoextra)
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
library(leaflet)
#Question 1
bb = read_csv('../data/uscities.csv') %>%
filter(city == "Palo") %>%
st_as_sf(coords = c("lng", "lat"), crs = 4326) %>%
st_transform(5070) %>%
st_buffer(5000) %>%
st_bbox() %>%
st_as_sfc() %>%
st_as_sf()
## Parsed with column specification:
## cols(
## city = col_character(),
## city_ascii = col_character(),
## state_id = col_character(),
## state_name = col_character(),
## county_fips = col_double(),
## county_name = col_character(),
## county_fips_all = col_character(),
## county_name_all = col_character(),
## lat = col_double(),
## lng = col_double(),
## population = col_double(),
## density = col_double(),
## source = col_character(),
## military = col_logical(),
## incorporated = col_logical(),
## timezone = col_character(),
## ranking = col_double(),
## zips = col_character(),
## id = col_double()
## )
#Question 2
bbwgs = bb %>% st_transform(4326)
bb = st_bbox((bbwgs))
q2_image = getlandsat::lsat_scenes() %>%
filter(min_lat <= bb$ymin, min_lon <= bb$xmin) %>%
filter(max_lat >= bb$ymax, max_lon >= bb$xmax) %>%
filter(as.Date(acquisitionDate) == as.Date("2016-09-26"))
## Parsed with column specification:
## cols(
## entityId = col_character(),
## acquisitionDate = col_datetime(format = ""),
## cloudCover = col_double(),
## processingLevel = col_character(),
## path = col_double(),
## row = col_double(),
## min_lat = col_double(),
## min_lon = col_double(),
## max_lat = col_double(),
## max_lon = col_double(),
## download_url = col_character()
## )
write.csv(q2_image, file = "../data/palo-flood-scene.csv")
md = read_csv("../data/palo-flood-scene.csv")
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
## X1 = col_double(),
## entityId = col_character(),
## acquisitionDate = col_datetime(format = ""),
## cloudCover = col_double(),
## processingLevel = col_character(),
## path = col_double(),
## row = col_double(),
## min_lat = col_double(),
## min_lon = col_double(),
## max_lat = col_double(),
## max_lon = col_double(),
## download_url = col_character()
## )
files = lsat_scene_files(md$download_url) %>%
filter(grepl(paste0("B", 1:6, ".TIF$", collapse= "|"), file)) %>%
arrange(file) %>%
pull(file)
st= sapply(files,lsat_image)
## File in cache
## File in cache
## File in cache
## File in cache
## File in cache
## File in cache
s=stack(st) %>% setNames(c(paste0("band", 1:6)))
cropper= bbwgs %>% st_transform(crs(s))
r=crop(s, cropper)
#Question 3 True
par(mfrow=c(1,2))
plotRGB(r, r=4 ,g=3 ,b=2, stretch= "lin")
plotRGB(r, r=4 ,g=3 ,b=2, stretch= "hist")
Infrared
plotRGB(r, r=5 ,g=4 ,b=3, stretch= "lin")
plotRGB(r, r=5 ,g=4 ,b=3, stretch= "hist")
False Color
plotRGB(r, r=5 ,g=6 ,b=4, stretch= "lin")
plotRGB(r, r=5 ,g=6 ,b=4, stretch= "hist")
My Choice
plotRGB(r, r=2 ,g=4 ,b=6, stretch= "lin")
plotRGB(r, r=2 ,g=4 ,b=6, stretch= "hist")
Stretching determines the representation of each value. Histogram stretching tends to be lighter and brighter, while linear stretching seems to be darker
5 Different Masks
palette= colorRampPalette(c("blue","white","red"))
ndvi=(r$band5 - r$band4) / (r$band5 + r$band4)
ndwi=(r$band3 - r$band5) / (r$band3 + r$band5)
mndwi=(r$band3 - r$band6) / (r$band3 + r$band6)
wri=(r$band3 + r$band4) / (r$band5 + r$band6)
swi= 1 / sqrt(r$band2 - r$band6)
## Warning in sqrt(getValues(x)): NaNs produced
Plotting
plot(ndvi,col=palette(256), legend=FALSE)
plot(ndwi,col=palette(256), legend=FALSE)
plot(mndwi,col=palette(256), legend=FALSE)
plot(wri,col=palette(256), legend=FALSE)
plot(swi,col=palette(256), legend=FALSE)
Water Thresholds
thresholding1= function(x){ifelse(x <= 0,1,NA)}
thresholding2= function(x){ifelse(x >= 0,1,NA)}
thresholding3= function(x){ifelse(x >= 0,1,NA)}
thresholding4= function(x){ifelse(x >= 1,1,NA)}
thresholding5= function(x){ifelse(x <= 5,1,NA)}
flood1= calc(ndvi,thresholding1)
flood2= calc(ndwi,thresholding2)
flood3= calc(mndwi,thresholding3)
flood4= calc(wri,thresholding4)
flood5= calc(swi,thresholding5)
plot(flood1, col= "blue", legend=FALSE)
plot(flood2, col= "blue", legend=FALSE)
plot(flood3, col= "blue", legend=FALSE)
plot(flood4, col= "blue", legend=FALSE)
plot(flood5, col= "blue", legend=FALSE)
Kmeans
set.seed(1)
r_extract = getValues(r)
r_nna = na.omit(r_extract)
kmeans_r = kmeans(r_extract, 12)
## Warning: Quick-TRANSfer stage steps exceeded maximum (= 5882000)
fviz_cluster(kmeans_r, geom="point", data = r_extract)
thresholding_s= function(x){ifelse(x >= 0,1,0)}
flood_s= calc(ndwi,thresholding_s)
kmeans_raster = flood_s
values(kmeans_raster) = kmeans_r$cluster
plot(kmeans_raster, col = viridis::viridis(12))
Finding the correct KMeans cluster that represents our water pixels
com_table = table(values(flood_s), values(kmeans_raster))
which.max(com_table[2,])
## 12
## 12
kmeans_raster[kmeans_raster != which.max(com_table[2,])] = 0
kmeans_raster[kmeans_raster != 0] =1
plot(kmeans_raster)
Changing remaining NAs to 0 for analysis
thresholdings1= function(x){ifelse(x <= 0,1,0)}
thresholdings3= function(x){ifelse(x >= 0,1,0)}
thresholdings4= function(x){ifelse(x >= 1,1,0)}
thresholdings5= function(x){ifelse(x <= 5,1,0)}
floods1= calc(ndvi,thresholdings1)
floods3= calc(mndwi,thresholdings3)
floods4= calc(wri,thresholdings4)
floods5= calc(swi,thresholdings5)
fs = stack(floods1, flood_s, floods3, floods4, floods5, kmeans_raster)
plot(fs)
sum_fs = fs %>% sum()
plot(sum_fs, col = blues9)
(cellStats(fs, sum) * res(fs)^2) /1e6
## layer.1 layer.2 layer.3 layer.4 layer.5 layer.6
## 5.9994 6.4908 10.7451 7.6221 13.6809 7.9371
How many layers classified the marked area in the video as a flood?
aoi = st_point(c(-91.78967, 42.06290)) %>%
st_sfc(crs = 4326) %>%
st_sf() %>%
st_transform(st_crs(fs))
raster::extract(sum_fs, aoi)
## [1] 3
There are 3 layers that classified this area as a flood; it is assumed that some rasters missed the marking of the area as a flood